function LEs = MLE(a, b, c, x0, y0, z0, N, transient, delta)
    eps_val = 1e-4;  % 更宽松的防止除零阈值

    X = [x0; y0; z0];
    Q = eye(3);
    sum_LE = zeros(3,1);
    
    for i = 1:N
        % 加强除零保护
        for k = 1:3
            if abs(X(k)) < eps_val
                X(k) = sign(X(k)) * eps_val;
            end
        end
        
        [x_next, y_next, z_next] = SIMM(X(1), X(2), X(3), a, b, c);
        X_next = [x_next; y_next; z_next];
        
        % 检测数值异常
        if any(isnan(X_next)) || any(isinf(X_next))
            % 直接跳过本步，状态不变
            X_next = X;
        end
        
        J = Jaco(X(1), X(2), X(3), a, b, c, delta);
        
        Z = J * Q;
        [Q, R] = qr(Z);
        
        % 防止对角元过小导致 log(0)
        r_diag = diag(R);
        r_diag(abs(r_diag) < 1e-12) = 1e-12;
        
        if i > transient
            sum_LE = sum_LE + log(abs(r_diag));
        end
        
        X = X_next;
    end
    
    LEs = sum_LE' / (N - transient);
end